Entropy maximization in the force network ensemble for granular solids 
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A long-standing issue in the area of granular media is the tail of the force distribution, in particular 
whether this is exponential, Gaussian, or even some other form. Here we resolve the issue for the case 
of the force network ensemble in two dimensions. We demonstrate that conservation of the total area 
of a reciprocal tiling, a direct consequence of local force balance, is crucial for predicting the local 
stress distribution. Maximizing entropy while conserving the tiling area and total pressure leads to 
a distribution of local pressures with a generically Gaussian tail that is in excellent agreement with 
numerics, both with and without friction and for two different contact networks. 

PACS numbers: 45.70.Cc, 05.40.a, 46.65.+g 



In a granular system the interactions among individ- 
ual grains are dissipative [l| , so that an undriven system 
eventually jams into a static, mechanically stable con- 
figuration [2[. There are typically many jammed states 
consistent with a given set of macroscopic constraints, 
e.g. fixed grain number and fixed pressure or volume. A 
granular packing can move from one jammed configura- 
tion to another only under externally imposed fluctua- 
tions such as shaking or tapping, because grains are too 
massive to rearrange thermally [l| . Despite this nonequi- 
librium nature, Edwards proposed that the methods of 
equilibrium statistical mechanics might describe many 
material properties of jammed media [3[. A successful 
statistical mechanics of jammed systems would represent 
an important theoretical handle on the physics of granu- 
lar media. It would, for example, permit the calculation 
of grain scale statistical properties, e.g. volumes [H or 
stresses [H, 0, 0, HI , from a small number of global con- 




FIG. 1: A force network (a) on the periodic frictionless trian- 
gular lattice. Edges represent contact forces; larger forces are 
thicker, (b) A wheel move. Arrows represent changes to the 
forces on each grain, (c) The reciprocal tiling corresponding 
to (a). Larger forces map to longer lines, (d) A move in the 
tiling. Shaded area is being exchanged among tiles. 



straints. 

In this Letter we derive an analytic expression for the 
probability density p(p) of pressures on individual grains 
in the bulk, which characterizes the strikingly hetero- 
geneous force networks observed in granular solids @. 
We work within the force network ensemble of Snoei- 
jer et al. [10] and confirm our results with highly accu- 
rate numerics. The force network ensemble comprises all 
balanced force configurations with a fixed global stress 
tensor (a) on a fixed contact network beyond marginal 
rigidity. The ensemble is both a minimal model for the 
statistics of local stresses and a basic test for any statisti- 
cal mechanical theory of stress states in granular systems. 

There is as yet no clear consensus on the form of the 
distribution of local stresses in granular media. Of par- 
ticular interest is the large-stress tail, which early exper- 
iments found to be exponential when measured on the 
boundary More recent measurements in the bulk 

[IH, along with numerics 0, find distributions that 
bend downward on a semilogarithmic plot, suggesting 
faster than exponential decay. A number of proposed 
theories exploit an analogy to the microcanonical ensem- 
ble to arrive at a Boltzmann-like exponential tail [1, 0, 0] . 
These theories should in principle apply to the force net- 
work ensemble, but numerical simulations tailor-made 
to accurately sample large contact forces unambiguously 
show a Gaussian tail in the force network ensemble in 
2D [3]. As the present statistical mechanics approaches 
fail to describe simple models like the force network en- 
semble, they must be missing an important ingredient. 
We argue that local force balance is absolutely crucial to 
describe the correct stress statistics. In particular, we 
show that the pressure distribution in the force network 
ensemble directly follows from entropy maximization, but 
only when it respects a conserved quantity overlooked in 
previous theories. This conserved quantity follows from 
force balance at the grain scale, and leads to excellent 
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agreement with numerics for both small and large forces. 

Force network ensemble. — Snoeijer's ensemble is com- 
posed of all "force networks", i.e. sets of noncohesive 
contact forces on a fixed granular contact network, for 
which all N grains are in static force and torque balance 
(e.g. Fig. []Ji). For packings with more than a critical 
number of contacts per grain z c , there exist many bal- 
anced force networks. z c = 4 (3) for frictionless (fric- 
tional) 2D packings of disks [10[. All force networks 
on a contact network with the same (a) and local force 
and torque balance can be sampled uniformly by a series 
of Monte Carlo moves, termed "wheel moves" 
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Fig. [lb gives an example. In the ensemble each force net- 
work has an equal a priori probability (a flat measure), 
in the spirit of the Edwards ensemble [3] . We illustrate 
our approach in the specific case of the periodic friction- 
less triangular lattice of circular grains before expanding 
to frictional packings and different contact networks. 

As by construction the global stress tensor (a) is fixed 
[l6j |. the extensive pressure in the system V is conserved, 
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V = ^ Pi = const . 



(i) 



The sum runs over all grains and pi = Tr<Ji/2 is the 
pressure on grain i. We restrict ourselves to isotropic 
states, (a) ~ 1, so that V fully characterizes (a). 

Every force network, regardless of its (dis)ordering, co- 
ordination number z, or friction coefficient /i, has a cor- 
responding reciprocal tiling. The systems in Fig. and 
lc are a real and reciprocal pair; each grain corresponds 
to a tile. Each face of the tile corresponds to one of the 
grain's contact forces. The face is oriented at a tt/2 ro- 
tation to the force /, and its length is proportional to 
|/|. Because the grain is in static force balance, the faces 
form a loop enclosing the tile [17[. By Newton's third 



law, the tiles fit together with no gaps. 

Specifying the boundary forces on a packing establishes 
the boundaries of its corresponding tiling, and hence the 
tiling's total area. Fixing (a) in a periodic system has the 
same effect. Rearrangements of bulk forces, i.e. the wheel 
moves, correspond to local exchanges of area among tiles. 
The total tiling area is unaltered by wheel moves, and 
therefore the area A is conserved. That is, 



N 

A = ai = const . 



(2) 



The sum runs over all tiles and is the area of tile i [18[ . 
The conservation of A can be seen explicitly for the fric- 
tionless triangular lattice in Fig.[T]i. Area conservation is 
a global constraint that results from imposing local force 
balance. It holds for arbitrary force balanced packings in 
2D with fixed (a) or boundary forces. It plays a crucial 
role in determining the statistics of local stresses. 




FIG. 2: (a) Fitted (dashed) and numerical (solid) average area 
of a tile with perimeter p.(b-d) Theoretical (dashed) and nu- 
merical (solid) pressure probability distributions for the fric- 
tionless triangular lattice with N — 1840. The additional 
numerical curves in (b) are from N — 460 and 115. 



We scale the grain diameter in the triangular lattice 
such that the pressure on a grain is the sum of its z = 6 
normal forces. In these units the perimeter of a tile is 
equal to the pressure on the corresponding grain, mak- 
ing the pressure p a convenient measure of local stress. 
Though the force distribution p(f) is more commonly 
studied, we expect that p(p) and p(f) have similar tails. 

Entropy maximization. — Armed with the insight that 
the force network ensemble involves two conserved quan- 
tities, V and A, we explore their implications for the 
statistics of local stresses. While previous work has in- 
corporated the conservation of V or its equivalent, the 
conservation of A has heretofore been overlooked. We 
will show that the conservation of A has a dramatic ef- 
fect on the force network statistics. 

We calculate the probability density p(p) by maximiz- 
ing entropy while conserving V and A. Each force net- 
work corresponds to a set of pressures {pi}, i = 1 . . . N. 
We define P(p)uj(p)dp as the probability of finding a pres- 
sure p in the interval + cu(p)dp), where oo(p) is the 
density of states for pressures. The entropy S is the log- 
arithm of the number of ways of constructing force net- 
works with pressures {pi} consistent with P(p). In the 
thermodynamic limit p 



(P(p) \nP{p))uo(p)dp. 



(3) 



The experimentally accessible probability density p(p) is 
related to P(p) by p(p)dp = P(p)uo(p)dp. 
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It is important to note that weighting all force net- 
works equally does not correspond to a flat measure on 
the pressures, i.e. u(p) ^ const. The contact forces {fi} 
on a grain, i = 1 . . . z, can be taken as coordinates of its 
state space. We demand that a grain explore only the 
regions of the space corresponding to force and torque 
balanced, noncohesive forces. We assume that, subject 
to these constraints and prior to imposing entropy max- 
imization, the grain is equally likely to be in any of its 
allowed states; this amounts to neglecting correlations 
with neighboring grains [l9[ • The result is a density of 
states that goes as u(p) oc p v . The value of v depends 
on the grain's coordination number and the friction co- 
efficient. For the frictionless triangular lattice, v = z — 3. 

The entropy is maximized subject to Eqs. (pQ) and ([2]), 
as well as normalization of p(p). This leads to 



1 = / p(p) dp , (p) = V/N = / pp(p) dp , 
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and (a) = A/N = / (a(p)) p(p) dp . (4) 
Jo 

(a(p)) = J a p(a\p) da is the average area of a tile with 
perimeter p; p(a\p) is the conditional probability a tile 
has area a given perimeter p. By the method of La- 
grange multipliers the entropy-maximizing density sub- 
ject to Eqs. (|4]) is 



p(p) = Z V exp(-ap-j(a(p))) . 



(5) 



Without the constraint on tiling area we would have 
7 = and an exponential tail: Incorporating local force 
balance by means of the area constraint has qualitatively 
changed the form of the distribution. The Lagrange 
multipliers Z, a, and 7 are determined by substituting 
Eq. ([5]) in Eqs. (|4]). For frictionless systems a scaling 
argument shows that (a(p)) is quadratic in the thermo- 
dynamic limit [20|]; we write (a(p)} = c(a)(p/(p)) 2 and 
determine the constant c from numerics. Thus the prob- 
ability density p(p) has a generically Gaussian tail, as 
was shown numerically for p(f) 

We employ umbrella sampling [14] on a periodic fric- 
tionless triangular lattice with N = 1840 to numeri- 
cally determine p{p). From the sampled (a(p)), shown in 
Fig. [2^i, we extract c « 0.89. Fig.[2j3-d contains the corre- 
sponding probability density of Eq. ([5j) and its numerical 
counterpart. Theory and numerics are in excellent agree- 
ment, even for p(p) as low as 10 ~ 8 . The slight discrep- 
ancies can be attributed to finite size effects and spatial 
pressure correlations: due to force balance, neighbors of 
large p grains are more likely to be at high p themselves. 
Thus large pressures are less entropically favorable than 
suggested by neglecting correlations. 

Frictional lattices. — We now consider frictional trian- 
gular [z = 6) and square [z = 4) lattices. A system 
with friction coefficient p permits contact forces with a 
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FIG. 3: (a,b) Theoretical (dashed) and numerical (solid) pres- 
sure probability distributions for 7-grain clusters in the fric- 
tional triangular lattice with p = 0.1, 0.5, 1.0, and 3.0. (c,d) 
Pressure distributions for 9-grain clusters in the frictional 
square lattice with p = 0.5, 1.0, and 2.0. (b, d inset) 7 of 
Eq. J5} for for various friction coefficients p. 



tangential component t < pn, n being the normal com- 
ponent. This has two important consequences. The first 
is that (a(p)) is not strictly quadratic. Friction permits 
tiles with area a < 0, which occurs when tile faces over- 
lap. Nevertheless, on dimensional grounds we expect 
(a(p)) ~ p 2 for large p. In numerics, deviations from a 
quadratic form increase with /i, but for all frictional sys- 
tems we have studied quadratic scaling holds for p > (p). 
Hence Eq. ([5]) still yields Gaussian tails. Secondly, we 
find that friction increases spatial correlations [l9[ • Con- 
sequently, as in Ref. [7[ , we coarse-grain and study clus- 
ters of m = 7 (9) grains and k = 30 (24) contacts on the 
triangular (square) lattice. The frictional clusters have 
exponent v = 2k — 3m — 1 in their density of states. We 
find (a(p)) for a cluster deviates much less from quadratic 
behavior than its single-grain counterpart. 

Lacking the exact form of (a(p)) for frictional systems, 
we determine the Lagrange multipliers satisfying Eqs. (jlj) 
using the numerically sampled (a(p)). Theory and nu- 
merics are again in excellent agreement, as seen in Fig.[3j 

Infinite friction. — As the Lagrange multiplier 7 tends 
towards zero with increasing p (Fig. [3)o and d, insets), 
we investigate the limit p — > 00. For finite friction and 
circular grains, normal and tangential forces are coupled 
through the force balance constraints on each grain and 
the Coulomb constraint on each force. In the infinite 
friction limit the Coulomb constraint has no effect. For 
the triangular lattice there are three distinct contacts 



per grain, and it is possible to choose a set of tangential 
forces {ti} to balance force and torque on each grain re- 
gardless of the normal forces {n^}. The only constraints 
on the {rii} are positivity, ni > 0, and fixed total pres- 
sure V. This leads directly to a Boltzmann distribution 
p(n) = exp (— n/(n)) with (n) = V/zN. In con- 

trast, for systems with z — z c < 3, such as the square 
lattice, the normal and tangential forces remain strongly 
coupled through force balance even for infinite friction. 
We have confirmed numerically that in the infinite fric- 
tion limit the Boltzmann distribution holds for the tri- 
angular lattice, and that the normal force and pressure 
distributions in the square lattice remain Gaussian. 

Boundary forces. — To this point we have imposed a 
flat measure on ensembles of periodic force networks. We 
have also numerically investigated the frictionless trian- 
gular lattice with a boundary and subject to a flat mea- 
sure on the boundary forces [21[. This produces by pre- 
scription a Boltzmann distribution of boundary forces, 
reminiscent of experiment [ll|. Nevertheless, we find 
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that the force and pressure distributions in the bulk, de- 
fined as grains at least six layers from the boundary, have 
Gaussian tails. This simple example demonstrates that 
a boundary distribution may not provide direct informa- 
tion about the bulk distribution. 

Conclusion. — We have derived an analytic expression 
for the distribution of pressures in the force network en- 
semble in 2D and found excellent agreement with numer- 
ics. Distinct from previous studies, we incorporate two 
conserved quantities, a total pressure and the area of a 
reciprocal tiling. The latter is a direct consequence of 
force balance on the grain scale, and we conclude that 
this is crucial to understand the statistics of local forces 
in granular media. As a result, large stresses obey Gaus- 
sian statistics. This observation is robust to changes in 
the contact network, the finite friction coefficient, and 
the imposed measure. 

We have not addressed the distribution at the unjam- 
ming transition, which could have a signature in the lo- 
cal stress statistics. Marginally rigid packings cannot be 
studied within the force network ensemble. Similarly, our 
results are restricted to two dimensions. A naive exten- 
sion of the reciprocal tiling to 3D suggests p (p) ~ e~ p 
with S = 3/2, while numerics find S « 1.7 [14| within 
the force network ensemble. The discrepancy may be the 
result of stronger spatial correlations than in 2D, where 
coarse-graining suffices, or it may signal new physics. 

Importantly, along with recent experiments [l, Il2j |. 
our results give serious cause to doubt that exponential 
statistics are a generic property of jammed granular mat- 
ter. At the very least, more work is needed to distinguish 
bulk and boundary phenomena and to clarify why mea- 
sured boundary forces show exponential statistics. 
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